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Using an inverse of the standard linear congruential random number generator, large randomly 
. occupied lattices can be visited by a random walker without having to determine the occupation 
[ status of every lattice site in advance. In seven dimensions, at the percolation threshold with L'' sites 
and L < 420, we confirm the expected time-dependence of the end-to-end distance (including the 
corrections to the asymptotic behavior). 
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, For transport in disordered media, like the ant-in-the-labyrith model of random percolation 

] [1, 2, 3], one usually first constructs the disordered medium, and then starts the transport 
^ ' medium. In a hypercubic lattice of linear size L in d dimensions, one first has to decide about 
the status of L'^ sites. In two dimensions this is efficient since for long enough time every site 
is visited. In higher dimensions, however, only a small fraction of the lattice is visited and 
^ it is more efficient to determine the status of a site only when it is visited first, keeping this 
^ . status fixed for later visits. We achieve this aim here for the case that every site is randomly 
Q ■ allowed for the random walker with probability p and forbidden with probability 1 — p. We 
apply this method to seven-dimensional percolation [4, 5] at its threshold Pc = 0.088951 where 
an infinite network of allowed sites is just possible [6]. We found good agreement of the new 
method [L < 420) with the traditional one {L = 23), while for smaller L finite-size effects are 
[ seen. 

Q ] Lattice sites i = 1,2, ... L are numbered in helical order such that the neighbours of site i 
^ ' are z±l, i ± L, i ± L'^, . . . ,i ± L'^~^. In the normal technique, we fill one site after the other 
using consecutive random numbers < Xi < 1; if Xi < p the site is allowed. With a simple 
^ . linear congruential random number generator like 

ji+i = kji mod m; m = 2" 

for integers ji instead of real numbers Xj, on computers with at least n bits per word, this 
method means that a step of the random walker to the right {i ^ i + 1) corresponds to a 
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multiplication of ji by k, a step to the left to a multiplication [7] by the inverse k* such that 
kk* = 1 modulo m. This inverse exists if and only if k and m axe relatively prime; then k* can 
easily be calculated by the extended Euclid algorithm [8]. Tables of * were given by 

L'Ecuyer [9] or can be easily calculated by a Fortran program available from DO; for example 
k = 16807 works with k* = 1278498327 for n = 32 and k* = 4409460005719528983 for n = 63. 
We used on a 64-bit Cray-T3E 

k = 3512401965023503517, k* = 3753721746144068021, 
table 5 of [9]. If a step upwards, i ^ i + L, corresponds to L multiplications with k (modulo m), 
then a downwards step corresponds to L multiplications with k* (modulo m), i.e. L pseudo- 
divisions by k. The other directions involve higher powers of L. In this way an arbitrary walk 
can be followed in a reproducible way by the proper number of multiplications with k and k*, 
always using the modulo restriction. The number of sites is then only limited by the period of 
the random number generator; we used a multiplier [9] with good spectral results and maximal 
period 2^^. This method can be extended to the more general generator jj+i = kji + c modm. 

We use helical boundary conditions in d — 1 directions and unlimited extent in the remaining 
direction, i.e. site i is thought to have as its right neighbour the site i + 1 even if i is an integer 
multiple of L and lies at the right boundary. And if due to a large jump the new site has 
i < or i > L"^ we do not have to put i back into the traditional interval 1 < i < L*^ through 
i ^ i±L'^ since the new method does not actually store an array with index i from 1 to L*^. In 
this sense our lattice size is infinite in one direction. (It is easiest to imagine a planar lattice, 
with sites numbered in a typewriter fashion.) 

In our 32-bit Fortran program, instead of the modulo function we used the automatic 
omission of leading bits if a product of two integers gives an integer with more than 32 bits. 
Since Fortran does not have the unsigned long integer type of C, the resulting products between 
—2^^ and +2^^ are often negative. We tested by comparison with the traditional method 
that the new technique works nevertheless. For 64 bits, we masked off the leading bit, thus 
working modulo 2^^ with only positive integers. This program including all shmem commands 
for communication between different Cray processors still fits onto 68 lines. Fig. 1. It can also 
be used for less than seven dimension by reducing the parameter idim; we made test runs with 
1,100,000,0002, l,050,000^ 32800^ 4100^ 1001^ sites, as well as 180^ 

c diffusion in d < 9 dimensions on random lattice; < 10"18 to avoid repeats 

implicit none 

real p,f actor ,av(30) ,avs(30) ,time 

integer L, idim,maxtime,num, i , iwalker ,maxt , iter , itime 
integer idir , isum, number , node 

integer*8 ishift(0: 15) ,mult(0: 15) ,ip,imod,iseed,isee2,ibm,ibmj ,new 
integer shmem_n_pes, shmem_my_pe, info, iadd, barrier 
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common /t3e/ av, avs 

data L/^2Q/ , idim/7/ ,p/0 . 088951/ ,maxtime/17/ ,num/1000/ 

data av/30*0 . 0/ , iseed/50923/ , isee2/471 1/ 

data imod/'7FFFFFFFFFFFFFFF'x/ 

node=shmem_my_pe 

number=shmem_n_pes () 

t ime=irtc 0*3. 33e-9 

if (node . eq. 0) print *,p,L,idim,iseed,maxtime,num, number 

ibm =2*iseed-l 
ibmj=2*isee2-l 
do 1 iter=0,node 
ibm =ibm *65539 

1 ibmj=ibmj*65539 

mult (0) =3512401965023503517 
mult ( 1 ) =375372 1746 14406802 1 
do 2 i=2,idim*2-l 

2 mult(i)=iand(mult(i-2)**L,imod) 
ip=p*imod 

f actor=l . 0/ (number*num) 

do 8 iwalker=l ,num 
do 3 i=l,idim 
ishift(2*i-2)=0 

3 ishift(2*i-l)=0 

4 ibmj=iand(ibmj*16807, imod) 
if (ibmj .gt . ip) goto 4 
maxt=2 

do 8 iter=l ,maxtime 

if (iter. gt. 2) maxt=maxt*2 
c evaluation of distances only if time is power of two 

do 6 itime=l,maxt 

5 ibm=ibm*16807 
idir=ishf t (ibm , -60) 

if (idir.ge.2*idim) goto 5 
new=iand(ibmj*mult(idir) ,imod) 
if (new.gt . ip) goto 6 
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ibmj=new 
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ishif t (idir) =ishif t (idir) +1 
continue 
isum=0 

do 7 i=l,idim 
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isum=isum+(ishift(2*i-2) -ishif t(2*i-l))**2 
av( iter) =av (iter) +isuin 



8 



continue 



info = barrier 
if (node.gt.O) stop 
if (number . gt . 1) then 
do 9 iadd = l,number-l 

call shmem_get (avs , av.mcLXtime , iadd) 
do 9 i=l ,ma:xtime 

9 av(i)=av(i)+avs(i) 
endif 

do 10 iter=l ,maxtime 

10 print *, iter, av (iter) *f actor ,2**iter 
time = irtcO * 3.33e-9 -time 

print *,'CPU= ' ,4*time, number 

stop 

end 



Figure 1: Fortran program for Cray-T3E parallel computing, averaging over number proces- 
sors each simulating num random walkers. Without the complications of parallel computation, 
loops 1 and 9 as well as all commands involving shmem, barrier, node, number, avs can 
be omitted. 

Now we repeat the standard scaling theory [1, 2, 3, 6] for the average end-to-end distance 
r{t,p) and enlarge it by predicting correction terms. 

Let Rs, < rs{t) >, < r{t) > be the average radius of gyration of a cluster containing s 
sites, the average end-to-end distance travelled on such a cluster in time t, and this distance 
averaged over all starting points on all clusters. Analogous notations hold for the squared 
distances; these averages < r^(i) > over the squares do in general not scale like the squares of 
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the averages < r(t) >. The number of clusters containing s sites each is = s~^/[(pc ~ p)s'^] 
with the standard critical exponents a = l/{f36), r = 2 + 1/5. In general dimensions d, we 
have Rs = s'^^h[{pc — p)s'^] with the correlation length exponent u, and assume 

< r',{t) > = R', g{t/s\ Rjs'^) . 
Right at p = Pc this assumption simplifies to 

<Tl{t)>^RlG{t/s^) , 

(varying as R^ for t ':$> s'^, when every site on the cluster has been visited many times, and as 
^2au/x _ similar to the anomalous behavior on the infinite cluster, for 1 ^ i ^ s^). For 

p < Pc and for sufficiently long times one has only the former behavior: the random walker has 
visited the whole finite cluster, and thus rs{t — oo) ~ Rg. Here, /, h, g, G, gi, g2, ■■■ are suitable 
scaling functions. 

Thus below Pc for long times we have 

< > ^^Rssnsgi[{pc-p)s''] 

s 

< r^{t) > ^^Rlsn,g2[{pc-p)s''] 

s 

in d dimensions. 

For d > 6 the critical exponents should be those of the Bethe lattice, with r = 5/2, a — 
V — 1/5 — 1/2, Us oc exp[— const (pc — pYs\. Also, = 6, x — 3/2. Then 

< r{t) > Const + . . . (pc - p)^^"^ 

< r'^{t) > oc log(pc — p) + . . . 

p < Pc- Right at p — Pc the powers oi Pc — p are replaced by the proper powers of 
-pf: 

<r{t)> ^ Const + 0(1 /t^/^) (1) 
< r^(t) > oclogt 

More precisely, in < r^(t) > oc s~^(y'(t/s^, const) we can approximate the sum by Z^^s"^ 
with an upper limit for s of order t^^^, giving ln(const t) or Const + In t. 

For the third-leading term in < r^(i) > at p = Pc, we note that the leading correction 
to scaling for d > 6 comes from the leading irrelevant parameter w, which represents the 
probability of having vertices with three bonds [10]. This implies corrections to any leading 
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power law behavior of the form (1 + const x w'^X^~'^), where X represents an appropriate length 
scale, which could be r(t), Rg or t^/'^^ . Specifically, at pc we expect corrections to < r^(t) > 
of relative order t^^~^^/^^ which become at c? = 7, just as for < r{t) > above. The two 

corrections differ for d > 7. Thus finally at (i = 7 we expect: 

< r'^{t) > oclogt + Const + 0{l/f/^) (2) 

Obviously, Eq. (1) is easier to test and reasonably confirmed by Fig. 2, using pc = 0.088951 
from [5]. For Eq. (2) we see in Fig. 3 a logarithmic long-time behaviour for L = 23 (traditional 
method) and L = 420 (new method), while for short times and/or small lattices deviations 
exist. Figure 4, similar to Fig. 2, shows that < r^(t) > — 1.25 ln(i) — 7.8 + S/t^^^ is consistent 
with our data for large systems over 8 decades in time. (< r'^{t — 1) — p exactly, but the curve 
shown in Fig. 4 extrapolates to a value ai t = 1 differing from this value by 0.2,) The fit is not 
shown in Fig. 3 since it would barely be distinguishable there from the curve for L = 420. 
Less extensive simulations in eight dimensions, L = 180, nicely confirm < > oc log(t) + 
const, but are not accurate enough to distinguish between correction exponents 1/3 and 1/6. 
In contrast, for seven dimensions 1/6 fits over a wider range than 1/3. 

In summary, the new method gives results consistent with the old one but allows for enor- 
mously larger lattice sizes; and these results are consistent with scaling theory. 

We thank Humboldt Foundation and GIF for supporting this collaboration, and R.M. Ziff, 
P. L'Ecuyer and P.M.C. de Oliveira for encouraging discussion. 
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Average distance versus correction-to-scaling power of reciprocal time, L = 7(+). 17(x), 20( 
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Figure 2: Average distance < r > versus l/time^^^; asymptotically a straight line is expected 
leading to a finite intercept. 
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L = 7 (', 300,000 walks), 17 (x. 10.000 walks'!. 23 64.000 walks) and 420 (curve. 64.000 walks) 
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Figure 3: < r^{t) > versus time on logarithmic time scale; L = 7 and 17 deviate for long times 
from L — 23 (traditional method) whereas L — 420 (continuous curve, new method) agrees 
with L — 23 for the observed times. 



<r r> -y*ln(t) versus correction-to-scaling power for L=420; y=1 .3, 1 .25, 1 .2 from top to bottom 
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Figure 4: < r'^{t) > -y\n{t) versus 1/time^/^ for y = 1.2, 1.25 and 1.3; (640,000 walks); the 
line comes from 64,000 walks to eight times longer t up to 128 Megasteps. Asymptotically a 
straight line is expected leading to a finite intercept. 
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